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We present a theory describing trapping of the normally dispersive radiation by the Raman solitons 
in optical fibers. Frequency of the radiation component is continuously blue shifting, while the soliton 
is red shifting. Underlying physics of the trapping effect is in the existence of the inertial gravity- 
like force acting on light in the accelerating frame of reference. We present analytical calculations 
of the rate of the opposing frequency shifts of the soliton and trapped radiation and find it to be 
greater than the rate of the red shift of the bare Raman soliton. Our findings are essential for 
understanding of the continuous shift of the high frequency edge of the supercontinuum spectra 
generated in photonic crystal fibers towards higher frequencies. 
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I. INTRODUCTION 

Frequency conversion in optical fibers has been an ac- 
tive research field already for few decades [1]. Most 
striking and extensively studied recent advance has been 
generation of extremely broad optical spectra (supercon- 
tinua) in optical fibers with small effective area, pumped 
by femto-second pulses with the carrier frequency close 
to the point of the zero group velocity dispersion (GVD) 
[2j [3]. Applications of the supercontinuum include 
spectroscopy, metrology [4 , telecommunication ^ and 
medicine [6]. 

Amongst problems posed by the observation of super- 
continuum one of the most puzzling has been under- 
standing of the nonlinear processes leading to the gen- 
eration of the high frequency wing of the supercontin- 
uum continuously drifting towards even higher frequen- 
cies [3 [HI m Uni [11] . Several experimental and numeri- 
cal observations explicitly demonstrated that the radia- 
tion at the blue wing of the supercontinuum often propa- 
gates in the form of nondispersive wave packets, localized 
on the femtosecond scale and continuously blue shifting 
Hmg. Note that GVD at the blue edge of the 
continuum is typically normal, therefore the dispersive 
spreading can not be compensated by the nonlinearity. 
Independently from the supercontinuum generation the 
effect of the localization of blue shifting pulses in the nor- 
mal GVD range, coupled to the Raman solitons propa- 
gating in the anomalous GVD range, has been reported 
in the series of papers by Nishizawa and Goto [161 [III [H] 
and more recently by Cheng and co-authors p^. It has 
been proposed in Refs. [151 Ull [HI that the physical 
mechanism behind the above effects is the cross-phase 
modulation (XPM) [1, 20 . 

However, it is well known that the XPM coupling be- 
tween anomalously and normally dispersing components 
can lead to dispersion compensation and formation of 
the bright-dark soliton pairs only if one soliton compo- 
nent is a dark pulse and the other one is bright, see, e.g., 
[21, 22, 23 . Thus the XPM can not be the sole reason for 
formation of the bright-bright localized states across the 
zero GVD wavelength [El [Tl [Ml [15l [H Ull [H] . Also, 



the red shift of the anomalously dispersing component 
is readily explained by the intrapulse Raman scattering 
[TJ [34], while the blue shift of the normally dispersing 
bright pulse coupled to it requires to be understood. 

Reference [TT has explained formation of the blue edge 
of supercontinua in fibers using the theory of four-wave 
mixing between the solitons and dispersive waves [24l [25] . 
It has been demonstrated that for typical fiber disper- 
sions the interaction between the soliton and the blue 
radiation happens recurrently, so that every scattering 
event leads to the further blue shift of the signal pulse 
[T4] . Though this theory describes well first stages of 
the blue edge formation it fails to explain why the fem- 
tosecond pulses emerging there remain free of dispersive 
spreading. The latter is naturally expected because of 
the strong normal GVD and would lead to a fast degra- 
dation of any nonlinear interaction between the soliton 
and the blue radiation, which in practice continues over 
long distances. In our recent work ^26^ we have explained 
the physical mechanisms behind existence of the non- 
dispersive and continuously blue shifting localized states 
of light on the high frequency edge of the supercontin- 
uum spectra. The light is trapped by the refractive in- 
dex changes induced, on one side (front of the pulse), by 
the red shifting Raman solitons via the nonlinear cross 
coupling and, on the other side (trailing tail), by the 
inertial force originating from the fact that the solitons 
move with acceleration. The nature of the latter effect is 
analogous to the gravity-like inertial force acting on an 
observer in a rocket moving with a constant acceleration. 

The aim of this work is not only to provide mathe- 
matical details for the mostly qualitative description pre- 
sented in [26 , but also to extend the theory into the 
regime of sufficiently strong intensities of the blue radia- 
tion. In this regime the trapped blue component of the 
two frequency bound-state starts to influence the soli- 
ton dynamics on the red edge, which makes noticeable 
quantitative impact on the propagation dynamics of the 
bound state. 
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FIG. 1: (Color online) Frequency dependence of the group 
velocity normalized to speed of light in vacuum (group index) 
(a) and of the GVD (b) in the photonic crystal fiber as in 
Ref. [14 . Fiber dispersion like this is typical for numerous 
supercontinuum experiments. 



II. THE MODEL 

The subject of this work is the explanation of the 
existence and detailed study of the previously reported 
[lailllllllSlIiailllllHlP^ two-frequency bright-bright 
soliton-like states in optical fibers, with frequency of one 
component being in the anomalous and frequency of the 
other being in the normal GVD range. Thus the disper- 
sion we need to use should include the sign change of the 
GVD. The simplest example of the dispersion operator 
having this property is 

d,E - \edlE = . (1) 

Assuming E — e*^^~*^^, we find that the GVD is = 
e^. Thus (5 measures the spectral deviation from the zero 
GVD point at ^ = 0. The third-order dispersion co- 
efficient e is positive for the telecom fibers and in the 
proximity of the 800nm for the typical photonic crystal 
fiber (PCF) designs used for generation of supercontin- 
uum with femtosecond pulses [2^ 3 . e can also be neg- 
ative in the proximity of 1500nm in PCFs and tapered 
fibers with sufficiently small cores [27]. It is important 
to note that the group index parameter d^h = e^^/2 is 
symmetric with respect to (5 —b and therefore the 
group velocity is matched across the zero GVD point. 
The matching, or near matching, of the group velocities 
is important for the existence of the two-frequency bound 
states, because the pulses should not spatially separate 
before the bound state is established. In the PCFs the 
third order dispersion d\k usually changes its value and 
sign as frequency varies, as can be seen from the change 
of slope of the GVD curve in Fig. [ijb). This eventually 
destroys the matching of the group velocities across the 
zero GVD point, see Fig.jlja). However, the fact that the 
matching is still achieved over the wide bandwidth is the 
most important for the effect of radiation trapping and 
for the formation of the blue wing of a supercontinuum. 

Taking into account the instantaneous Kerr nonlinear- 
ity and the non-instantaneous Raman response, the light 



propagation in a fiber is modeled by the dimensionless 
generalized nonlinear Schrodinger (NLS) equation pLj 

d^E = ik{idt)E + i{l - 0)E\Ef 

/ + 00 
Rit')\E{t-t',z)fdt'. (2) 
-CXO 

The dispersion operator in Eq. ([2| is 

^ r Rirn) 

Hidt)^Y.^i'dr, (3) 

m=3 

where the coefficients /S'^^^ are selected to fit the disper- 
sion profile of the fibers. The total electric field is given 
by ^e*^r2;-za;r-t -^q^q^^ where the reference frequency cOr is 
chosen to coincide with the zero GVD frequency, that is 
why the sum in Eq. (|3| starts from m = 3. R{t) is the 
standard Raman response function of silica: 

R{t) = e ^^^l^ Q(t)e-'/^^ sin - . (4) 

T1T2 Ti 

Here 9(t) is the Heaviside function and parameter = 
0.18 weights the Raman nonlinearity relative to the Kerr 
one. Characteristic times of the delayed Raman response 
are ri = 12.2/s/r and r2 = 32/s/r [1 . t is the dimen- 
sionless time in the reference frame moving with the light 
group velocity at coq and measured in the units of r. z 
is the distance along the fiber measured in the units of 
L, where L is any convenient characteristic length. Field 
amplitude E is scaled to where 7 is the nonlin- 

ear parameter of the fiber [1 . To get feel for the real 
values of the parameters we choose 7 = 0.02 (Wm)~^ 
and L = 8m, which gives 1/(7^) = 6.25W. Choosing 
r = 200fs and = 0.12ps^ /km [corresponding to the 
dispersion slope at the zero GVD point close to 800nm 
in Fig.[l];b)] we have e = Lp^^^r^ = 0.12. 

III. NUMERICAL EXPERIMENTS 
ILLUSTRATING RADIATION TRAPPING BY 
ACCELERATING SOLITONS 

We proceed describing two sets of numerical experi- 
ments. First, is when supercontinuum evolves from a 
single pump pulse. Second, is when the two-frequency 
bound state is excited directly by the two pulses. In the 
first set we used the realistic fiber dispersion, see Fig. [l] 
However, to explain the effect of radiation trapping by 
solitons across the zero GVD point it is sufficient to con- 
sider the simple cubic dispersion as in Eq. ([T]), which 
simplifies comparison of analytical and numerical results 
and is used throughout the rest of the paper. 

Figs. [2] illustrate supercontinuum generation with a 
single pump in the fiber as in Fig. [l] The pump fre- 
quency was chosen either in the range of anomalous 
GVD, Figs. [2|a),(b), or in the range of normal GVD 
Figs. [2](c),(d). In both cases one can see that, after the 
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FIG. 2: (Color online) Left panel: spectral evolution along 
the fiber length for the dispersion profile as in Fig.[l] Pump 
pulse parameters: (a) wavelength 900nm is in the anomalous 
GVD range, (c) 800nm is in the normal GVD range. In both 
cases peak power is 6kW, duration 200fs; The radiation at 
the short wavelength edge of the continuum starts its contin- 
uous drift towards shorter frequencies at z > 0.2m for (a) and 
z > Im for (c). Right panel: spectrograms showing simulta- 
neous frequency and time domain images of the supercontinua 
in (a) and (c) for z = 1.5m and z = 3m, respectively. The 
function plotted has been calculated using the XFROG inte- 
gral /(t,cj)= J^^dt'Eref{t'-t)E{t')e-''^'\ ReieEref is the 
2ps gaussian pulse. See Ref. [28] on more details about the 
XFROG (cross-correlation frequency resolved optical gating). 



initial stages of the supercontinuum development (see e.g. 
[14 , for detailed description), the blue tip of the spec- 
trum starts its continuous drift towards higher frequen- 
cies, see z > 0.2m in Fig. |2|a) and z > Im in Fig. [2|c). 
This spectral shift appears to be correlated with the soli- 
ton self-frequency shift at the opposite (infrared) edge of 
the spectrum. Spectrograms showing the signals simulta- 
neously in the frequency and time domains, see Figs.[2|b) 
and (d), unambiguously demonstrate that the high fre- 
quency tip of the continuum is localized in the time do- 
main on the same femtosecond scale as the soliton at the 
infrared edge. One can also see, that not only the soliton 
at the very edge of the spectrum, but all the red shift- 
ing solitons have associated localized pulses on the high 
frequency side of the spectrum in the normal GVD range. 

To isolate the effect of formation of the two-frequency 
bound-states across the zero GVD point we perform sim- 
ulations using the two pulse excitation. From now on 
we will focus on the simple case of cubic dispersion as in 
Eq. ([T]). We start from the case e > 0, see Fig.js] The first 
of the pump pulses is spectrally located in the anomalous 
GVD regime and forms a red shifting Raman soliton. The 
second pulse is delayed with respect to the first one and 
has spectrum in the range of normal GVD. In the course 
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FIG. 3: (Color online) Spectral (a) and time-domain (b) 
evolution along the fiber pumped with two pulses obtained 
from the numerical modehng of Eq. ^ with p'^'^>^'> = 
and e = +0.12. Initial conditions are E(z = 0,t) = 
^2^1 sech( t/n) exp(— 'i^it) + A/2^sech((t— To)/t2) exp(— 2^2^). 
with qi = 2000 (peak power 24:kW), q2 = 250 (peak 
power ^ 3kW), ri — 

y/-D'J/2q2 ^ 0.57 ( 
Si -80 (- 964nm), S2 80 (- 684nm). Results 
in (b) are presented in the acceleratin g fr ame of reference: 
t' — t — ^o^^/2, with qq defined by Eq. ( 15 ). 



D'^/2qi = 0.2 (- 40/5), T2 = 
115/s), To = 0.25 (- 50/s) and 



of the propagation the latter pulse appears to be trapped 
on the trailing tail of the soliton, while its frequency is 
continuously increasing, see Fig. [3| Note that the Ra- 
man effect pulls the soliton toward smaller frequencies 
and away from the zero GVD point, therefore the GVD 
felt by the soliton is continuously increasing. This leads 
to the noticeable temporal broadening of the soliton and 
to the drop in its amplitude, which also affects the shape 
of the trapped radiation, see Fig.[3];b). Clearly the group 
velocities of the two components within the bound state 
are the same, which suggests that the frequencies are 
changing along the dispersion curve in such a way that 
the group velocity matching is preserved. 

In the case of the negative third-order dispersion (e < 
0) the effect of the radiation trapping is also clearly 
observed, see Fig. [4j The difference here is that the 
normally and anomalously dispersing components of the 
two-frequency bound-state are now spectrally converging 
towards the zero GVD point, see Fig.|4ja). However, the 
direction of the spectral evolution of each of the compo- 
nents is the same as for e > 0. Namely the Raman effect 
pulls the soliton component in the anomalous GVD range 
towards smaller frequencies, while the pulse in the nor- 
mal GVD range shifts towards higher frequencies. The 
GVD felt by the soliton is reducing in this case, therefore 
the soliton is adiabatically compressed, see Figs. |4]^b), 
(f). When the soliton is pulled sufficiently close to the 
zero GVD point, so that the significant part of its own 
spectrum is in the normal GVD regime, it starts to emit 
Cherenkov radiation [29l|30l[3l], see Figs.|4|;a), (e). This 
radiation creates spectral recoil effect, which counter- 
balances the Raman self-frequency shift and stabilizes 
the soliton frequency [2TI |32] . As the soliton looses its 
energy through radiation, it broadens and so does the 
trapped wavepacket. Still, both spectral components of 
the bright-bright quasi-soliton pair remain fairly well lo- 
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FIG. 4: (Color online) Pulse propagation in the case of neg- 
ative third-order dispersion, e = —0.12. (a) and (b): two 
pulses pump similar to that in Fig. [3] Parameters of the 
pulses are: q\ — 200 (peak power =~ 2AkW), q2 — 25 
(peak power - O.SkW), n = -D'l j^qx ^ 0.63 (- 130/s), 

T2 = y^-D'J/2q2 ^ 1.8 (- 360/s), To = 1.5 (- 300/s) and 
Si = 80, 62 = -80. (c) and (d): XFROG spectrograms for 
the case of two-pulse pump, calculated at different propaga- 
tion distances, (e) and (f): soliton propagation with q = 200 
(peak power ~ 2AkW). All time domain results are presented 



in the accelerating frame of reference: t' — t 
go defined by Eq. (15). 



goz /2, with 



calized over propagation distances much larger than the 
GVD length. Spectrograms showing the effect of the ra- 
diation trapping in this case are shown in Figs, [i] (c), 
(d). Note, that the stabilization of the frequency of the 
pure soliton, i.e. without the trapped radiation, happens 
at longer propagation distances, cf. Figs. [4]^a) and (e). 
This suggests that the trapping effect boosts the rate of 
the soliton self frequency shift. We will discuss this more 
in details in Sec. VI. 



IV. COUPLED NLS EQUATIONS IN THE 
ACCELERATED FRAME OF REFERENCE 



In order to explain the above numerical observations 
of the two-frequency bound states across the zero GVD 
point, we reduce Eq. ([2| to the coupled NLS equations 
for the two pulses on the opposite sides of the zero GVD 



frequency. We assume 

E = Ai ex.-p[ikiz — i5it] + A2 exp[z/c2^ — ^^2^], (5) 

where Ai^2 are the amplitudes of the two pulses, 5i^2 are 
their frequencies and ki^2 = ^(^1,2) are the wavenum- 
bers: k{S) = e^^^k{idt)e~'^^^ . We also assume that 81^2 
are selected in such a way that k'^ = k2 {k' = d^k) and 
therefore the group velocities of the two pulses are equal. 
Substituting Eq. (5| into Eq. ([2|, expanding \E{t-t')\^ 
in Eq. ([2| up to the first order Taylor term and neglecting 
all the fast oscillating exponential terms with frequencies 
^1^2 and their harmonics we obtain a pair of the coupled 
NLS equations: 



id^Ai + didlAi 
id^A2 + d2dlA2 



\A2\'] 



-[IA1P + 2IA2P] A, 
TAid, [\A, 
-[IA2P + 2IA1P] A2 
TA2d4\A,\'^\A2\'] 



(6) 



(7) 



where x = {t-zk[)/ y^\k^\,T = is the effective 

Raman time, Td = J tR{t)dt ^ 0.0073, di = 1/2 (anoma- 
lous GVD), and ^2 = -\k'^\/\2k'{\ < (normal GVD). 
The most important restriction of Eqs. ([6|, ^ relative 
to Eq. ([2| is that the former do not include frequency 
dependence of the GVD, see Fig. [6) However, this is not 
critical for understanding of the trapping mechanisms. 

If A2 = then Eq. (|6| has an approximate solution in 
the form of the NLS soliton moving with constant accel- 
eration g [1 . So that its center in the (x, z)-plane follows 
the parabolic trajectory x = gz'^/2 and its frequency is 
continuously red shifting with the rate g/{2di). For the 
single NLS equation it has been demonstrated that there 
exist a transformation into the accelerating frame of ref- 
erence, which retains the structure of the NLS equation 
apart from adding a linear in x potential [33 . ^ is a free 
parameter of this symmetry transformation, which value 
is fixed by assuming that the corrections to the soliton 
induced by the linear potential are balanced by the cor- 
rections due to the Raman term, see below. 

Here we apply the analogous symmetry transformation 
to the coupled equations 



Ai = V^(2;,0exp 
^2=0 {z, exp 



^2d^ +^/u^) 



(8) 
(9) 



where 



^ = x- xo{z), xo = gz^/2, /i,2(^) = --r^ i ^o^^- 

4ai,2 J 

(10) 

Parameters q and A are the shifts of the propagation con- 
stants of the two components. 

The resulting equations for and (j) are 

id.ij^d^dl^ = #- [|V^|' + 2|(/)|2] + 



5 



id. 



T0S^ [IV^P + 1^1'] + • (12) 



9^ 

2d2 



The prime difference of the coupled NLS equations (11), 
(12) with the text book ones [1 is the presence of the 



hnear in ^ potentials. The acceleration g has been the 
free parameter upto now. Its selection will be discussed 
in the next two sections. 

We should note here, that the natural next step in 
our analysis could be the setting up a boundary value 
problem for the independent version of Eqs. ( 11 ), ( 12 ). 



However, proper setting of the boundary conditions is a 
challenging problem due to presence of the linear poten- 
tial. The main reason for this is that working out asymp- 
totic behavior of solutions for ^ ^ — oo requires a delicate 
analysis going beyond the scope of this work. Indeed, 
if one simply neglects the nonlinearity, then the tails 
of both components behave like Airy functions, which 
was assumed in the prior works on the similar problems 
[35l[36l[^. However, the amplitude of the oscillatory tail 
— oo) of the Airy function decays only as |^|~^/^, 
which makes some linear terms in Eqs. (11), (12) to de- 



cay at the rates matching the decay rate of the nonlinear 
terms, suggesting that equating the solution tails to the 
Airy function is an approximation. This problem still 
waits for its proper analysis even in the case of the single 
soliton equation. For the above reasons we rely in what 
follows on perturbation theory complemented by the nu- 
merical calculations with zero boundary conditions for 
large |^|, which correctly describes the localized parts of 
the solutions, but gives only qualitative answers at the 
oscillatory tails. The effects we study below are, however, 
dependent most strongly on the localized part of the so- 
lutions, and therefore our approximation is adequate in 
the given context. 



V. LINEAR THEORY OF THE RADIATION 
TRAPPING 

We proceed, by considering the limit when the nor- 
mally dispersive component (j) is much weaker than the 
soliton component ip. In this limit we can assume that 
the field t/^ is not affected by and neglect all the terms 
nonlinear in Then Eq. (11) for ip can be solved per- 
turbatively. We assume = ?/^o(0 + V^ilO + • • - ^ where 



2(7sech(V'g/di^) 



(13) 



solves did?ipo 



t^u — qijjQ — ipQ and ipi accounts for corrections 
due to Raman effect and the linear potential: 



didjipi - #1 



3V^o'V^i=TV^oS^V^o' + ^V^o. (14) 



The operator did'^ — <? + 3?/^o self-adjoint and singular. 
Its null space is spanned by the single eigenfunction d^ijjQ. 



Projecting the right-hand side of Eq. (14) on the latter 
we find 



9 = 90 



32Tq'^ 
15 



(15) 



Thus the linear potential indeed compensates for the Ra- 
man shift, at least in the first order, providing g = go- 
Note, that this approach neglects the non-exponential de- 
cay of the soliton tail at ^ ^ — oo, see discussion at the 
end of the previous section, and it is equivalent to the 
traditional considerations 031 [34] , where the oscillatory 
tail has been ignored. 

The ^^-independent equation for <p is simply a linear 
eigenvalue problem in our approximation 



\d2\dl<p + v{0<P = x<l>, 



(16) 



where potential V{^) consists of the localized soliton part 
and of the linear potential induced by the acceleration 
with already determined value: 



ViO = 2i>l - Td^i>l + 



2M2I 



(17) 



The superposition of the exponentially decaying soliton 
tail and of the linear potential creates a local minimum of 
V on the trailing tail of the soliton, see Fig. [sja), giving 
rise to the effect of light localization in the normal GVD 
regime. The soliton itself serves as a potential barrier for 
the normally dispersive waves on one side of the well and 
the linear potential, obtained as a result of the transfor- 
mation to the accelerated frame of reference, serves as a 
barrier on the other side. Thus the accelerating potential 
creates an inertial force acting on photons. It is analo- 
gous to the gravity-like force acting on massive bodies in 
a closed container moving with a constant acceleration. 
The soliton created potential wall is not infinite, how- 
ever. Therefore the light is expected to tunnel through 
it. Hence, only quasi-bound states embedded inside the 
continuum are possible as solutions of Eq. (16). 



In order to locate these states we first replace the soli- 
ton part of V by its asymptotic and find that 



V 
A 



(18) 



16^, B = 2 



C 



90 

2\d2 



Contribution to A due to the Raman term (~ T) is small 
as compared to IQq and is neglected. The potential Vh 
goes to infinity on both sides and has the discrete set 



of true bound states with eigenvalues An\ Eq. (|16[) 



with V replaced by Vh has been solved numerically and 
some of its eigenstates are shown in Fig. [5] with dashed 
lines. Then we take the finite potential V and for each 
n find few eigenvalues and corresponding eigenstates in 
the spectral proximity of \n \ applying the zero bound- 
ary conditions at both ends. The latter implies that we 
reliably calculate only the states with relatively small tail 
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FIG. 5: (Color online) (a) Full (dashed) line shows the poten- 
tial V (Vb) for the soliton with q = 100 and T = 0.0024. (b-e) 
show selected quasi-bound eigenstates of the potential V, cf. 
horizontal levels in (a). Dashed lines show the corresponding 
modes of the potential Vb- 



amplitude at ^ — oo. Having attempting more precise 
calculations would go beyond our original level of preci- 
sion anyway, because the oscillations of the soliton tail at 
^ ^ — oo have been disregarded in the first place. Some 
quasi-bound eigenstates of the true potential V are shown 
in Fig.|5|b)-(e) with the full lines, and the corresponding 
eigenvalues are indicated in Fig. [sja). 

The terms exp [^i;o/(2(ii,2)] in Eqs. (|8|, ^ explicitly 
express the continuous frequency shifts of the two com- 
ponents with the rates g/2di^2' The anomalous GVD 
{di > 0) corresponds to the expected red frequency shift 
[34] . However, the normal GVD {d2 < 0) implies the blue 
frequency shift. The latter explains spectral dynamics of 
the trapped radiation. Note, that some analogy may ex- 
ist here with the previously reported blue shift of the 
spectral holes associated with the dark fiber solitons ex- 
isting in the normal GVD range [38]. 

Equations ([8|, ([9| give, however, little physical in- 
sight into which elementary wave scattering mechanisms 
lead to these opposing frequency shifts. The physical 




FIG. 6: (Color online) Numerical propagation of the first 
three bound states within coupled equations (|6|, ([7|. Initial 
conditions are: soliton for Ai with ^ = 100 and linear eigen- 
state of V for A2. Left column: spectrum of the total field E, 
Eq. (Isl, with S2 — —5i — 20. Right column: temporal profiles 



of \Ai\'^ (dashed lines) and \A2\ (solid lines). Time domain 
results are presented in the acceleratin g fr ame of reference: 
t' — t — qqz^ with qq defined by Eq. ( 15 ). 



process driving the red shift of the soliton component 
is the well known intrapulse Raman scattering [TJ [3^ . 
The blue shift of the radiation component is driven by 
the intrapulse four-wave mixing described in details in 
[14 . Briefly, it means that the scattering of the radiation 
pulse on the soliton generates the blue shifted pulse. The 
continuous frequency shift of the soliton and the phase 
matching conditions work out in such a way that the ob- 
servable result of this process is the continuous blue shift 
of the radiation pulse carrier frequency. Trapping effect 
sustains this process over the long propagation distances 
and results in existence of the stationary soliton-radiation 
states. 

In order to verify validity of our approximate eigen- 
value analysis we have initialized Eqs. (|6|, ^ with the 
soliton for Ai and the linear eigenstates of V for A2 and 
solved the equations numerically. Fig. [6] shows evolution 
of the first three quasi-bound eigenstates. The amplitude 
of the eigenstates has been kept small in order to ensure 
that we remain in the regime, when 1^2^ terms are neg- 
ligible. One can see that, despite using the simplified 
boundary conditions, our solutions satisfy well the cou- 
pled NLS equations. This simulation also confirms that 
the propagation distances on which the tunneling induced 
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losses lead to noticeable effects are much larger than the 
typical GVD length. So that for the fiber length of or- 
der meters the dispersive spreading of the high frequency 
radiation is suppressed and it propagates as a localized 
state of light. Note also, that the characteristic feature 
of each quasi-bound state, apart from the lowest one, is 
the presence of several spectral peaks. Such multi-peak 
spectral structures are typical for the blue wing of super- 
continua seen in Figs.[2|a), (c). 

Approximation V ^ Vb is also useful because it allows 
us to carry out explicit variational calculations of the 
eigenvalues and eigenfunctions, and thus to have analyt- 
ical estimates for the width of the trapped states. Let 
us consider the variational approximation for the ground 
state (n = 0) only. As a trial function we choose 



6o = exp[-(C-7)V"''] 



(19) 



where 7 is the shift of the intensity maximum of the 
trapped state with respect to the soliton one and w is 
the width of trapped state. One can suggest a better 
trial functions accounting for the asymmetry of the pro- 
file of the ground. However, our choice is the best suited 
for getting transparent analytical expressions for param- 
eters 7 and w. The variational estimate Aq for the true 
value Ao < Aq is 



-\d2\dl 



(20) 



Minimizing Aq with respect to 7 and w'^ we find 

= 

7 = 
Ao = 




(21) 
(22) 
(23) 



Thus the narrow solitons (large q) and large Raman ef- 
fect (large T) result, quite naturally, in stronger local- 
ization of the radiation. Figure [7|^a) illustrates compar- 
ison between the numerically calculated linear mode 0o 



of the potential V{^) (T7|) and the variational approxi- 
mation. The numerical] 



ly calculated dependencies Ao(^) 
and Aq^^((7) for potentials V and 14, respectively, compare 
well with the variational result (23), see Fig. [7|^b). 



VI. NONLINEAR THEORY OF THE 
TWO-FREQUENCY QUASI-SOLITONS ACROSS 
THE ZERO GVD POINT 

The theory presented above explains the nature of the 
radiation trapping and is adequate in the regime when 
the radiation is weak. In this regime the family of the 




100 



FIG. 7: (Color online) (a) Lowest trapped linear mode of 
the potential V (solid line) and it s va riational approxima- 



tion (dashed line), given by Eq. (19) with 7 and w de 



fined through (22) and (21), respectively. Soliton parameter 
is ^ = 100. (b) Aq^'*(^), solid line, Ao(^), circles, and Ao(^), 
dashed line. 



soliton-radiation bound states is continuously parame- 
terized by q only, i.e. by the amplitude of the soliton 
pulse. In particular, this is expressed in the fact that the 
acceleration parameter g = go and the eigenvalue A are 
fixed for a given q. In the nonlinear regime one should 
expect that A will become a continuously varying param- 
eter, as it happens for other types of localized solutions 
in incoherently coupled NLS equations [27. Thus the 
acceleration g is expected to be continuously parameter- 
ized by both q and A, i.e. by the energies of the both 
fields (j) and t/j. In order to demonstrate explicitly, that 
indeed there is a problem to be addressed here, we com- 
pare the propagation distance at which the spectral recoil 
from the Cherenkov radiation stabilizing the soliton fre- 
quency takes place with and without the blue detuned 
pulse seeded into the fiber. One can see. Figs. [4]^a) and 
(e), that the sharp transition to the regime without the 
self-frequency shift happens at a shorter distance with 
the blue radiation present. This is because the spectral 
peak of the soliton reaches the critical distance from the 
zero GVD point sooner, which indicates that the rate of 
the self-frequency shift is faster for the soliton-radiation 
bound state than for the pure soliton, i.e. g > go- Thus 
nonlinear in (f) corrections should be taken into account 
to explain this effect. 

One obvious small parameter in our problem is the Ra- 
man time T, which enters both equations for (j) and ip. 
The linear in ^ potential can be considered as a small per- 
turbation only in the equation for i/j. However, it plays 
a crucial role in the localization of the ^ component at 
^ ^ 00 and therefore can not be neglected already in the 
leading order in the equation for (j). As we have found 
above, see Figs. |5] and [6) the (j) and components overlap 
only by their tails. Therefore the nonlinear coupling can 
be considered as a small perturbation on the compo- 
nent, which is self- localized. However, in the (j) equation 
the term is the only localization mechanism for 

^ ^ — 00 and hence can not be neglected there. Now we 
rewrite Eqs. (11), (12) collecting all the leading terms in 



the left hand-side, all the first order corrections in the 
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right-hand side and neglecting the rest: 



2di 



-2\4>\^xP + Txpd^\ij(' 



id,^-\d2\dl^-\<p+\<p\^<p + 2\i>\^^ 
= T<^d^\4>\'' + ... 



2d2 



(24) 



(25) 



In order to calculate the deviation of g from qq^ we as- 
sume that ^ = V^o(0+^i(0 + - ••A = </>o(0 + </>i(0 + - • 
where ipi and 4>i have the same order of smallness as the 



right-hand sides in Eqs. ( |24[ ), ([25|. This leads to 
where (po obeys 



2dx 



^0 - '^fyo, 



(26) 



(27) 



Projecting the right-hand side of Eq. (26) on d^ijjQ. 
We find 




FIG. 8: (Color online) Numerically calcul ated acceleration g 
for different branches of solutions of Eqs. (27), (28) with the 



potential V being replaced by its asymptotic Vb. Circles cor- 
respond to the solutions in the full potential V . Approximate 
analytical solution (30) is indicated by dashed line. Insets 
illustrate solution profiles (asymptotic potential) at different 
values of the parameter A. 



d 

— / d(^(t)l^l)Qd^i)Q, 
Q J —00 



(28) 



One should remember though that 0o itself is a function 
of g. Therefore Eq. (28) is an equation for g^ which 



needs to be solved. Solving Eq. (27) numerically we find 
the set of functions 0o parameterized by g and A, which 
are the direct continuations of the linear discrete set of 
eigenfunctions of the potential V found in the previous 
section. The tail of the solution for ^ ^ —00 is oscillatory 
and weakly decaying like the tail of the linear solutions. 
However, despite the fact that ^0 is only semi-bound, 
it can be used in the integral ([28|, because it is multi- 
plied there by the exponentially localized functions 
and d^tpO' It refiects the fact that the tail of ^ has only 
a negligible contribution into the selection of g. There- 
fore, like in the previous section, we can replace V with 
V5 and carry out calculations using the infinite potential 
well and exponentially decaying solutions. Substituting 
00 inside the condition (28) and solving the latter nu- 
merically for ^, we find the corrected values of g. Fig. [S] 
shows dependencies of ^ on A for the first three bound 
states. One can see that increase in the amplitude of the 
0-component generally leads to the larger values of ^, 
which practically means stronger negative accelerations 
and larger frequency shifts (the red/blue shift for the 
i/j / ^-component , respectively) . 

It is useful to derive an approximate analytical expres- 
sion for ^, which can be done in several ways. First, 
the variational approach can be applied to the nonlinear 
problem ( [27| ). However, this leads to the rather cum- 
bersome and difficult to understand expressions. More 
elegant answer, which also matches well our numerical 



calculations in Fig. [Sj can be obtained in the limit of the 
weak nonlinearity in Eq. (27). This is accomplished by 
solving Eq. (27) perturbatively. We consider the ground 
state and assume that in the ffrst order (j)o and A are given 
by their variational approximations Acf)^ and Aq found in 
the linear case, see Eqs. (19), (20). Here A is the con- 



stant amplitude to be determined. We also introduce cor- 
rection to the eigenfunction, 5(j)^ and to the eigenvalue, 
(5A = A — Ao, induced by the nonlinearity. Taking e as a 
dummy small parameter, we assume A'^ ~ e, Sc/) ^ e and 
SX ~ e. Looking at the integral in Eq. (28) we see that 
it has order e^, because its value is proportional to A'^ 
and it acquires an extra order of smallness due to a small 
overlap between 0o and i/jq. It means that g = go-\-0{e'^). 
The resulting equation for 5(j) derived from Eq. (27) is 



(■ 



-\d2\dl + V{^,go) 



\o)S<p ^ 



5XA4>o 



(29) 



The operator in the left hand side is self-adjoint and sin- 
gular. Its null- space is given by the linear ground state. 
Therefore projecting the right-hand side on (j)o we ffnd 
SX = V2A'^. Substituting 0o in Eq. (28) with A^o and 
taking into account that inside the integral in Eq. (28) 



the ipo can be approximately replaced with its asymp- 
totic 2^y2qexp[—^^yq/dl] we ffnd a relatively simple ex- 
pression for g: 



g = go 



,5/8 



A — Aq d\ 



151^2 
32 



1/4 



(30) 

One can see that the above assumption about smallness 
of the overlap between (po and ^Jjq remains valid provided 
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FIG. 9: (Color online) Numerical propagation of the low- 
est bound state within coupled equations ( |11[ ) and (12) with 
acceleration g = qq^ (a) and (b), and g = g, (c) and (d). 
Left column: anomalous GVD component (Ai), right col- 
umn: normal GVD component {A2). Intensities are plotted 
in logarithmic scale, (e) and (f): maximum of the intensity 
of Ai and A2 components, respectively, along propagation 
distance. Initial shape for A2 was calculated with the ap- 
proximate potential Vb for the soliton parameter ^ = 100 and 
A = 70 (^A ^ 20). Ai was initialized with the soliton (13) 
with q = 105 in order to account for initial radiative losses, 
see text for more details. 



^-9/82^-1/4 ^ ^ ^-^j^ 2^1/4 ]^g-^g practically order of 1, 
it implies that q should be sufficiently large. Together 
with the variational approximation for Aq, see Eq. (23), 



the above equation agrees very well with numerically cal- 
culated dependence of ^ on A for the ground states, see 
Fig-i 

An ultimate method to confirm the validity of our ap- 
proximate calculations is to take Eqs. (11), (12) in the 



reference frame moving with acceleration g = go and with 
g = g given by Eq. (30), initialize them with the bound 



state and solve numerically. Then trajectories of the so- 
lutions on the (z, ^)-plane should be straight lines in the 
case of the exactly selected acceleration and parabolic 
otherwise. Fig. |9] demonstrates the results of this numer- 
ical experiment, where intensities of the two components 
Ai and A2 are plotted separately. Since initial excita- 
tion was taken in a localized form, i.e. without proper 
account of oscillating tails, one can see that both com- 
ponents emit radiation during propagation. The radia- 



tion leakage from the trapped component, A2, is small 
[note logarithmic scale in Figs. [9] (a)-(d)] and does not 
considerably affect the acceleration of the bound state. 
More important, however, is the initial outburst of ra- 
diation from the soliton component, Ai. As the result, 
a noticeable drop in intensity of Ai (~ 5%) is observed 
during initial stage of soliton propagation, see Fig. [9](e), 
while the intensity of the A2 component stays practically 
the same during propagation, see Fig. [9](f). In order to 
account for these losses, in the initial condition for Ai 
we have increased the soliton parameter q hy 5% with 
respect to the value used in the calculations of the accel- 
eration g. Taking this into account, the Eq. (30) gives 



a very good approximation for the acceleration of the 
bound state and, hence, for the rate of the self-frequency 
shift associated with it, see Figs, ^c) and (d). 



VII. SUMMARY 

In summary, we have presented the detailed theory of 
the effect of radiation trapping by the Raman acceler- 
ated fiber solitons responsible for formation of the blue 
wing of the supercontinuum spectra. We demonstrated, 
that the radiation in the range of the normal GVD is 
subject to the inertial force, which, together with the 
soliton induced refractive index change, forms an effec- 
tive potential well prohibiting the dispersive spreading 
of the radiation. We have found not only the ground 
state of the radiation field, but also its excited states, 
which relevance for the past and ongoing experimental 
observations is under current investigation. The soliton- 
radiation bound states move with a constant acceleration 
in the time-space, while in the spectral domain the peaks 
corresponding to the two components move in opposite 
spectral directions (soliton component always gets red- 
der, while the radiation component gets bluer). 

In the first part of our theoretical considerations we 
have assumed that the radiation is linear. In this case 
we have demonstrated that the continuous blue shift of 
its frequency happens at the same rate as has been pre- 
viously calculated for the red shifting solitons J] [34] . 
Considering effects nonlinear in the radiation, we have 
found that the acceleration of the soliton-radiation bound 
states increases with the radiation amplitude, which cor- 
responds to larger rates of the self- frequency shift. We 
have also derived an approximate analytical expression 
for the latter. 

The results presented above pave the way for design of 
new fiber based soliton frequency converters, which al- 
low for efficient blue frequency shifts of the soliton like 
state. This removes traditional restriction of the Raman 
soliton based frequency conversion been directed only to- 
wards longer wavelengths. Our results also emphasize 
that moving dielectric media can be created via nonlin- 
ear modulation of the refractive index by the pulse prop- 
agating with a speed of light. This creates an interesting 
testing bed for studies of the effects of light propagation 
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in moving dielectrics, which have generated significant 
recent interest, see e.g. f39]. 
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